{--
    It was observed that the 3 year numbers 2013, 2014 and 2015 
    each have 3 prime factors that are pairwise different.
    
    The question was raised when this will happen next time.
    
    Definitions: 
    
    - A natural number n > 1 is _k-kinded_, when it is the product of 
    exactly k prime factors and every factor occurs 
    exactly once in the product.
    
    - A triple (a,b,c) is _interesting_ if b is the successor of a and
    c is the successor of b and a, b and c are k-kinded for some k.
    
    - A _triple of rank k_ is an interesting triple where the members are k-kinded.
    
    Computations:
    
    - A list of all interesting triples.
    - A list of all triples of rank k, for some given k
    - The smallest triple of rank k for some k
    
     Observations:
     
     - Exactly one member of an interesting triple will be a multiple of 3.
     This is because in 3 successive numbers, exactly one is divisible by 3.
     
     - The first number of an interesting triple is of the form 4n+1. This is
     self-evident as soon as one realizes that no number divisible by 4 can be
     k-kinded for any k, since it contains the factor 2 twice. Hence, an interesting
     triple can only occupy the 3 numbers between 4n and 4(n+1), for some n.
     
     Open Questions:
     
     Does a triple of rank k exist for every k>1?
     Are there infinitely many triples of rank k for some k?
     Are there infinitely many triples of rank k for all k?
     
     We know for sure that triples of rank 1 do not exist, as prime number triples
     are impossible. Empirical evidence seems to show that triples of rank 2 are
     quite abundant, while higher ranked triples appear much less often.
     
     Up to now, the smallest triple of rank 5 has been found at 16467033.      
    -}
    
module examples.YearNumbers 
            inline ( kprime, remi, rem9i, ° , °°, rem4s, rem9s, nextRem, nextRem1, nextRem3)
        where

import examples.EulerLib

type Kind = Int
type Rank = Int
type Zahl = Integer
type Faktoren = [Zahl]
type Zerlegung = (Zahl, [Zahl])
type Tripel = (Zerlegung, Zerlegung, Zerlegung)

--
--  Items defined elsewhere:
--
--  Zahl.primes     the list of prime numbers of type Zahl
--  length xs       the length of a list
--  null xs         true if and only if the list is empty, otherwise false
--  a `rem` b       remainder, optimized for natural numbers
--  a `quot` b      quotient, optimized for natural numbers
--  a:as            the list with head a and tail as
--  []              the empty list
--  [a]             shorthand for a:[]
--  [a,b]           shorthand for a:b:[]
--  succ a          successor of a
--  iterate f a     a : iterate f (f a)
--  filter p xs     all elements of list xs that satisfy predicate p
--  map f (a:as)    (f a) : map f as
--  (f . g) a       f (g a)
--  head (a:as)     a


--- Check if a number is k-kinded
kinded :: Zahl -> Kind -> Bool
kinded z k = g z k false Zahl.primes
    where
        -- Parameters:
        -- z    the number to check
        -- k    the kind we want
        -- u    is true, if the top prime was already a factor
        --      hence, if it divides z again, the result is false
        --      as the prime factors are not pairwise unique.
        -- p    the next prime number to check
        -- ps   subsequent prime numbers
        g :: Zahl -> Kind -> Bool -> [Zahl] -> Bool
        g !z !k !u (!p:ps)
            -- | traceLn("g " 
            --     ++ show z ++ "  "
            --     ++ show k ++ "  "
            --     ++ show u ++ "  "
            --     ++ show p)      = undefined 
            | k < 1             = false                     -- too many kinds
            | pot ==  z         = if k==1 then not u else false  -- z is itself prime or higher kinded
            | k == 1, p*p > z   = true              -- this is a prime             --               
            | k > 1, pot > z    = false             -- z has less than k factors
            | z `rem` p == zero = if u then false   -- p was already factor
                                    else g (z `quot` p) (k-1) true (p:ps)
            | otherwise         = g z k false ps
            where pot = p `hoch` k
        g _ _ _ _ = false  
        -- fast exponentiation function
        hoch :: Zahl -> Int -> Zahl
        p `hoch` k = h one p k
            where
                h :: Zahl -> Zahl -> Int -> Zahl
                h !a !p 0 = a
                h !a !p k = h (a*p) p (k-1) 

--- Check if a number is the head of a tuple of rank k
ranked :: Zahl -> Kind -> Bool
z `ranked` k = (z `kinded` k)  
                && (sz `kinded` k) 
                && (ssz `kinded` k)
    where
        sz  = succ z
        ssz = succ sz

--- The number 4
four :: Zahl
four = fromInt 4

--- The lower bound of the 'f' function
--- The smallest prime that consists of k odd prime factor
lowerbound k = (prod . take k . tail) Zahl.primes

--- Give the next number of the form 4n+1
succ4fold :: Zahl -> Zahl
succ4fold n = case n `rem` four of
    r | r == zero       = succ n
      | r == one        = n
      | r == succ one   = (succ . succ . succ) n
      | otherwise       = (succ . succ) n
      
--- The list of the first members of all k-ranked triples 
triples k =  [ m |  m <- [n, n+4 ..], m `ranked` k ]
    where !n = (succ4fold . lowerbound) k

--- Candidates for k-ranked triples
candidates k = (filter (`kinded` k) . iterate (four+) . succ4fold . lowerbound) k

--- The first member of the first tuple with rank k
f = head . triples 

--- Auxiliary code to print a nice result
result n = do
    let nfs n = show n ++ "=" ++ joined "*" (map show (factors n))
    println   (nfs n ++ ", " ++ nfs (succ n) ++ ", " ++ nfs (succ (succ n)))
                                                                              

badk _ = println "give a small number for k"
badn _ = println "give a number for n"

goodk k = do
    println ("lower bound is " ++ show (lowerbound k))
    println ("starting search at " ++ (show . succ4fold . lowerbound) k)
    let r = f k
    result r
    println ("f " ++ show k ++ " = " ++ show r)

goodn k n = do
    println ("lower bound is " ++ show (lowerbound k))
    println ("starting search at " ++ (show . succ4fold . lowerbound) k)
    mapM_ result (take n (triples k))

goodk2 sn k = either badn (goodn k) (String.int sn)

main ["kinded", sz, sk] = println ( fromInt sz.atoi `kinded` sk.atoi)
main ["walk", sz, sq] = do
        let k = sz.atoi
            q = sq.atoi
        println ("k=" ++ show k)
        println ("q=" ++ show q ++ ", p_q=" ++ show (kprime parr q)) 
        walk q ((foldr Liste Ende . reverse . take k . iterate succ) 0)

main ["candidates", sk, sz] = println 
    (length (takeWhile (< fromInt sz.atoi) (candidates sk.atoi))) 
main ["find", sk] = findK sk.atoi    
main ["check"] = mapM_ (cacheCheck (Cache.new 48)) [0..47]
main [sk] = either badk goodk sk.int
main [sk, sn] = either badk (goodk2 sn) sk.int
main _ = mapM_ println idee


idee = [ (p,q,r) |
    p <- Int.primes,
    p != 3, p != 11,
    b = p*33+1,             -- second element of tuple
    b `rem` 34 == 0,        -- divisible by 2*17
    q = b `quot` 34,
    q != 2, q != 17,
    q.isPrime,
    c = b+1,
    c `rem` 35 == 0,        -- divisible by 5*7,
    r = c `quot` 35,
    r != 5, r != 7,
    r.isPrime
        ]

cacheCheck :: Cache -> Int -> IO ()
cacheCheck cache i = do
    print ("index=" ++ show i ++ ", ")
    print ("prime=" ++ show (kprime cache.primes i) ++ ", ")
    print ("remainder=" ++ show (remi cache i) ++ ", ")
    println ""

{--
    Alternative idea.
    
    Instead of going through the numbers and checking whether 
    they are k-kinded, we build k-kinded numbers and check
    whether they are the head of a triple with rank k.
    
    To do this, we use a data structure called a _triangle_,
    that contains all the combinations of (k-1) prime numbers up to,
    but not including, p.
    
    For example, for k=3 and p=11, the triangle looks like this:
    
    > 11,3,2    11,5,2  11,7,2      -- products increasing to the right in each row
    >           11,5,3  11,7,3      -- products increasing downwards in each col
    >                   11,7,5
    
    (of course, for higher k, we get more dimensions.)
    
    The essential idea is that we do not need to actually store this triangle.
    Rather we just know how to traverse one. Note that this does not yield
    products in sequential order, though.
    
    In addition, we can easily construct the lowest and the highest
    member of each p-triangle. The lowest value is located in the top left,
    while the highest one is in the bottom right.
    
    To make things efficient, we do not store actual prime numbers but
    zero based indexes of prime numbers. The indexes are isomorph to the
    prime numbers with regard to relational operators.
    
    In addition, we do not store p, we just know the p of the triangle
    we are dealing with.
    
    Last, but not least, we completely ignore the first row in the search
    for candidates, as the elements there yield only even numbers.
      
    -}

data KI = Ende 
        | Liste {!top :: Int, !rest :: KI} -- (k-1) indexes in decreasing order, the first being the
                        -- one that is varied most often.
                        -- Invariant: The indexes further down in the list
                        -- are smaller than all previous ones.

type KP = [Zahl]   -- can be obtained from KI by replacing indexes with primes                                                  

parr = arrayFromList (take 1000 Zahl.primes)

infixl 6 `°` `°°`

(°) :: Int -> Int -> Int 
a ° b = if a == b then 1 else 3

(°°) :: Int -> Int -> Int
a °° b = (a*b) `rem` 9

--- lookup up prime number corresponding to i (O(1))
kprime :: JArray Zahl -> Int -> Zahl
kprime arr i = elemAt arr i

mapKI :: (Int -> a) -> KI -> [a]
mapKI f Ende = []
mapKI f Liste{top, rest} = f top : mapKI f rest

--- map a index list to a list of primes
kprimes :: JArray Zahl -> KI -> KP
kprimes  = mapKI . kprime

--- get the product of the primes corresponding to indexes, assuming we are in triangle q
kproduct arr q = reduce arr (kprime arr q) 
    where
        reduce !arr !p Ende             = p
        reduce !arr !p Liste{top, rest} = reduce arr (p * kprime arr top) rest
        


               
next :: Int -> KI -> KI
next !qi Ende = Ende
next !qi Liste{top, rest}
    | top+1 < qi = Liste (top+1) rest
    | otherwise = case next (qi-1) rest of
        Ende  = Ende
        liste@Liste{top=t} = Liste (t+1) liste 

--- remainder modulo 4 of a prime
remp cache p = (p `rem` four).int

--- remainder modulo 4 of the prime at index i
remi Cache{rems} i = rems.[i]

--- remainder modulo 4 of the prime at index i
rem9i Cache{rem9} i = rem9.[i]


--- remainder modulo 4 of a KI in triangle q
rem4s cache q = reduce4 cache (remi cache q)

reduce4 !cache !r1 Ende = r1 
reduce4 !cache !r1 Liste{top, rest} = reduce4 cache (r1°remi cache top) rest

--- remainder modulo 9 of a KI in triangle q
rem9s cache q = reduce9 cache (rem9i cache q)

reduce9 !cache !r1 Ende = r1 
reduce9 !cache !r1 Liste{top, rest} = reduce9 cache (r1°°rem9i cache top) rest

--- next index that has the same remainder as i (or 0, if outside cache)
nextRem Cache{same} i = same.[i]

--- next index that has remainder 3, or 0 if there is none
nextRem3 Cache{rem3} i = rem3.[i]    

--- next index that has remainder 1, or 0 if there is none
nextRem1 Cache{rem1} i = rem1.[i]

--- given a KI with remainder 1, compute next KI with remainder 1
nextR :: Cache -> Int -> KI -> KI
nextR !cache !q Liste{top, rest} 
    = case nextRem cache top of
        j | j > 0, j < q = Liste j rest       -- easy case
          | otherwise = normR cache q (next (q-1) rest)
nextR !cache !q Ende = Ende

--- Not so easy case - we must recompute the remainders.
--- Takes a (k-1) index list and returns a k index list!
--- Can also be used to normalize a tuple.
normR !cache !q Ende = Ende
normR !cache !q is = case rem4s cache q is of
    1 -> case nextRem1 cache is.top of
            j | j > 0, j < q = Liste j is       -- found one in next row
              | otherwise    = normR cache q (next (q-1) is)  -- try again
    3 -> case nextRem3 cache is.top of
            j | j > 0, j < q = Liste j is
              | otherwise    = normR cache q (next (q-1) is)  -- try again
    r -> error ("normR: how can remainder be " ++ show r ++ "?")

--- Go down to beginning of next row
downR :: Cache -> Int -> KI -> KI
downR !cache !q Liste{top, rest}  = normR cache q (next (q-1) rest)
downR !cache _ Ende  = Ende

--- Construct the KI that starts at line i of triangle q
--- i.e, for k=5 construct i+3,i+2,i+1,i
--- All products right/down from this one are greater. 
lineKI :: Int -> Int -> Int -> KI
lineKI k q i
    | i+k-2 < q = (foldr Liste Ende . take (k-1) . (iterate pred)) (i+k-2)
    | otherwise = Ende                    

--- Last line in triangle q
--- q-1,q-2,q-3,q-4
lastLine k q = q-(k-1)

--- Given a tuple, tell us the line
line Liste{top, rest}
    | Ende <- rest = top
    | otherwise    = line rest
line Ende = error ("line Ende")

--- Given a triangle, tell the last line number where numbers < p can be found
range !cache !k !p !q = r cache k p q 1 (lastLine k q)
    where
        r :: Cache -> Int -> Zahl -> Int -> Int -> Int -> Int
        r !cache !k !p !q !first !last 
            | midp > p     = r cache k p q first mid
            | mid+1 < last = r cache k p q mid last
            | otherwise    = if kproduct cache.primes q (lineKI k q last) > p 
                                then last-1 
                                else last
            where
                mid = (first + last) `quot` 2
                midp = kproduct cache.primes q (lineKI k q mid)     

--- Given a number and a triangle        
                        
walk qi t = do
    println (kprimes parr t)
    case t of
        Ende ->   return ()
        _    ->   walk qi (next qi t)

--- give us the lowest k-triple (excluding q, as usual)  
minki k q 
    | k <= q = (foldr Liste Ende . take (k-1) . iterate pred) (k-1)
    | otherwise = error "k > q?"

--- because 'minki' does not depend on q, all k-triangles have the same minimum
start k = minki k k
 
--- the greatest k-triple in the qtriangle (excluding q, as usual)
maxki k q
    | k <= q = (foldr Liste Ende . take (k-1) . iterate pred) (q-1)
    | otherwise = error "k > q?"                       

nine  = fromInt 9 :: Zahl
seven = fromInt 7 :: Zahl

--- Find first triple by sequentially walking a triangle with nextR
--- The argument must be an element whose product has remainder 1 modulo 4
findhitR :: Cache -> Int -> Int -> KI -> Int -> (Int, KI)
findhitR !cache !k !q Ende !fc = (fc, Ende)
findhitR !cache !k !q ki !fc 
    | rem9s cache q ki >= 7 = findhitR cache k q (nextR cache q ki) fc
    -- | traceLn("candidate " ++ show z) = undefined
    | z <- kproduct cache.primes q ki,
      succ z `kinded` k, succ (succ z) `kinded` k = (fc, ki) 
    | otherwise = findhitR cache k q (nextR cache q ki) (fc+1)
    

--- Find hit below p
findhitLT :: Cache -> Int -> Zahl -> Int -> Int -> KI -> Int -> (Int, KI)
findhitLT !cache !k !p !q !ln Ende !fc = (fc, Ende)
findhitLT !cache !k !p !q !ln ki   !fc
    | line ki > ln              = (fc, Ende)
    | z > p                     = findhitLT cache k p q ln (downR cache q ki) fc
    | rem9s cache q ki >= 7     = findhitLT cache k p q ln (nextR cache q ki) fc
    | succ z `kinded` k, 
      succ (succ z) `kinded` k  = (fc, ki)
    | otherwise                 = findhitLT cache k p q ln (nextR cache q ki) (fc+1)
    where
        z = kproduct cache.primes q ki

-- findhit :: Cache -> Int -> Int -> KI -> Int -> (Int, KI)
-- findhit !cache !k !q [] !fc = (fc, ende)
-- findhit !cache !k !q ki !fc 
--     | z `rem` four != one = findhit cache k q (next q ki) fc
--     | succ z `kinded` k, succ (succ z) `kinded` k = (fc, ki)
--     | otherwise = findhit cache k q (next q ki) (fc+1)
--     where
--         arr = cache.primes
--         z = kproduct arr q ki


{--
    We maintain a cache with 4 arrays of the same length:
    
    - array of prime numbers, 
    - the remainder modulo 4 for the prime number with same index
    - the index of the next prime number with the same remainder
      as the one at the index
    - the index of the next prime number with remainder 1
    - the index of the next prime number with remainder 3 
    
    The cache will get expanded as the need arises. We give initially 1000
    elements, and when a cache with larger size is requested, we compute one
    with the requested size rounded up to the next multiple of 1000, 
    to keep frequent cache recomputations in check.
    
     OTOH, we want the cache to be as small as possible, so that it'll fit in the
     CPU cache. 
    -}
data Cache = Cache {
                !primes :: JArray Zahl,
                !rems   :: JArray Int,
                !rem9   :: JArray Int,
                !same   :: JArray Int,
                !rem1   :: JArray Int,
                !rem3   :: JArray Int } where
    --- ensure cache stores element with given index 
    ensure :: Cache -> Int -> Cache
    ensure cache i
        | i < cache.primes.length = cache
        | otherwise = Cache.new (((i+1025) `quot` 1024)*1024)
    --- create a new Cache of required size
    new :: Int -> Cache
    new upper 
        | traceLn("expanding cache to size " ++ show upper) = undefined
        | otherwise = Cache { primes = ps, rems, rem9, same, rem1, rem3 }
        where
            ps = arrayFromList (take upper Zahl.primes)
            rems = (arrayFromList . take upper . map (Zahl.int . (`rem` four))) Zahl.primes
            rem9 = (arrayFromList . take upper . map (Zahl.int . (`rem` nine))) Zahl.primes
            findrem r ab
                | ab < rems.length = if rems.[ab] == r then ab else findrem r (ab+1)
                | otherwise = 0
            samerem i = findrem rems.[i] (i+1)
            same = arrayFromList (map samerem [0..upper-1])
            rem1 = arrayFromList (map (findrem 1) [1..upper])
            rem3 = arrayFromList (map (findrem 3) [1..upper])                     

{- 
mit Restklassen
findK k=8, q=36, min=761425665, max=126459568506372769
3153160534097993=157*149*127*109*107*101*53*17, 3153160534097994=238339*337*193*167*29*7*3*2, 3153160534097995=6121*2111*47*43*41*31*19*5
2562723 candidates checked without luck.
runtime 385.864 wallclock seconds.

ohne Restklassen
findK k=8, q=36, min=761425665, max=126459568506372769
3153160534097993=157*149*127*109*107*101*53*17, 3153160534097994=238339*337*193*167*29*7*3*2, 3153160534097995=6121*2111*47*43*41*31*19*5
2562723 candidates checked without luck.
runtime 459.002 wallclock seconds.

div Optimierungen, inkl. Restklassen 9:
findK k=8, q=36, min=761425665, max=126459568506372769
3153160534097993=157*149*127*109*107*101*53*17, 3153160534097994=238339*337*193*167*29*7*3*2, 3153160534097995=6121*2111*47*43*41*31*19*5
1932698 candidates checked without luck.
runtime 303.26 wallclock seconds.

-}
findK k = go (Cache.new 1024) k k
    where
        go :: Cache -> Int -> Int -> IO ()
        go cache k q = do
            let t0 = normR cache q ((start k).rest)
                -- t1 = start k
            println ("findK k=" ++ show k
                ++ ", q=" ++ show q
                ++ ", min=" ++ show (kproduct cache.primes q t0)
                ++ ", max=" ++ show (kproduct cache.primes q (maxki k q)))
            case findhitR cache k q t0 0 of
                (cs, Ende) -> do
                    println ("nothing found, " ++ show cs ++ " candidates checked without luck.")
                    go (cache.ensure (succ q)) k (succ q)
                (cs, t) -> do
                    let p' = kproduct cache.primes q t
                    result p'                
                    println (show cs ++ " candidates checked without luck.")
                    p <- down cache k p' q (nextR cache q t)
                    result p
        down :: Cache -> Int -> Zahl -> Int -> KI -> IO Zahl
        down cache k p q Ende = down cache k p (q+1) (normR cache (q+1) ((start k).rest))
        down cache k p q ki
            | min > p = return p
            | otherwise = do
                let maxln = range cache k p q
                println("findP k=" ++ show k
                    ++ ", p=" ++ show p
                    ++ ", q=" ++ show q
                    ++ ", ln=" ++ show maxln
                    ++ ", min=" ++ show min)
                case findhitLT cache k p q maxln ki 0 of
                    (cs, Ende) -> do
                        println ("nothing found, " ++ show cs ++ " candidates checked without luck.")
                        let cache' = cache.ensure (succ q)
                            t1 = normR cache' (q+1) ((start k).rest)
                        down cache' k p (succ q) t1
                    (cs, t) -> do
                        let p' = kproduct cache.primes q t
                        result p'
                        println (show cs ++ " candidates checked without luck.")
                        down cache k p' q (nextR cache q t) 
            where
                t0 = normR cache q     ((start k).rest)
                min = kproduct cache.primes q (start k)
